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Abstract 

Direct numerical simulation of open channel flow over a geometrically rough wall has been 
performed at a bulk Reynolds number of Ret ~ 2900. The wall consisted of a layer of spheres 
in a square arrangement. Two cases have been considered. In the first case the spheres are 
small (with diameter equivalent to 10.7 wall units) and the limit of the hydraulically smooth 
flow regime is approached. In the second case the spheres are more than three times larger 
(49.3 wall units) and the flow is in the transitionally rough flow regime. Special emphasis is 
given on the characterisation of the force and torque acting on a particle due to the turbulent 
flow. It is found that in both cases the mean drag, lift and spanwise torque are to a large 
extent produced at the top region of the particle surface. The intensity of the particle force 
fluctuations is significantly larger in the large-sphere case, while the trend differs for the 
fluctuations of the individual components of the torque. A simplified model is used to show 
that the torque fluctuations might be explained by the spheres acting as a filter with respect 
to the size of the flow scales which can effectively generate torque fluctuations. Fluctuations of 
both force and torque are found to exhibit strongly non-Gaussian probability density functions 
with particularly long tails, an effect which is more pronounced in the small-sphere case. Some 
implications of the present results for sediment erosion are briefly discussed. 

1 Introduction 

Sediment erosion by turbulent open-channel flow is an important aspect for fluvial engineering 
applications as it can, for example, cause the collapse of bridges. The mechanisms that lead to 
sediment erosion however are far from being completely understood due to the complex interactions 
between the turbulent flow and the sediment particles. Turbulent flow in an open-channel is 
statistically inhomogeneous in at least one spatial direction, and the Reynolds numbers of interest 
are typically high, which leads to a wide range of velocity and length scales. In addition, the 
presence of a range of sediment sizes, shapes and compositions complicates the description. Under 
certain conditions, the turbulent motions might erode the bed and entrain sediments as a result of 
the hydrodynamic force and torque acting on the particles. In order to improve the understanding 
of the mechanisms that lead to sediment erosion it appears necessary to simplify the problem 
under consideration. Therefore, as a first step we study the statistical properties of hydrodynamic 
force and torque acting on fixed spherical particles adjacent to the wall-plane in fully-developed, 
open-channel flow. For this configuration, we have performed direct numerical simulations (DNS). 

A wide body of literature exists that focuses on the hydrodynamic force acting on spherical 
objects placed in a fluid flow. In the low Reynolds number range analytical solutions have been 
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proposed for vario us flow configurations, e.g. the cas e of a particle in a linear shear flow (jSaffmanl . 
1965 ; Auton , 1987 ) , in a non- uniform rot ational flow ( Auton. Hunt fc Prud'homme , 19881) , and of a 



particle in the vicinity of a smooth wall (jKrishnan fe Leighton . 19951 ). In order to gain information 
on the mechanism that leads to lift and drag on a particle in the range from small to moderate 
Reynold s numbers, simila r flow configurations have also been explored by mean s of experimental 



studies (IKing fe Leightonl. 119971) and by means of direct numerical simula tions (|Kim et all [1993 



iBagchi fc Balachandarl 120021 : IZeng et all 120091 : iLee fc Balachandarl lioioh . In the high Reynolds 
number l imit numerous studies ca n be found that describe the flow around spheres in unbounded 
flow (see lYun. Kim fc Choi . 



2006, for an overview) 



The studies mentioned above have focused on situations in which the flow field approaching 
the sphere is laminar in nature. It is well known, however, that turbulence can have a significant 
effect on the statistics of the force s acting on a sphere. A revie w on the effect of turbulence on 
an isolated sphere can be found in IBagchi fc Balachandarl (|2003l ) . The authors studied the forces 
on an isolated sphere subject to free-stream isotropic turbulence for small and moderate Reynolds 
numbers by means of direct numerical simulation. They found that turbulence had only little 
effect on the mean drag and that the fluctuations of lift and drag scaled linearly with both the 
mean drag and the turbulence intensity. 

In contrast, turbulence appears to have a si gnificant effect in the case of a sphere positioned 
close to a wall, the lift being particularly affected ( Willetts fc Murravlll98ll : [Halll . 1988; Zcn g et al. 



20081 ). The experimental evidence shows, that similar to the low Reynolds number regime, sig- 
nificant positive values fo r mean lift (directed away from the wall) are obtained for a s phere 
touching the wall plane (IWilletts fc Murravl Il98ll lHalll . Il988t iMollinger fc Nieuwstadtl . Il996 : 



Muthanna. Nieuwstadt fc Huntl . 120051 ). When the sphere is not touching the wall, the picture 
is less cle ar and still a matter of d iscussion: both positive and negative values of the lift are 
reported. IWilletts fc Murray J 1981 ) fo und changes in sign for the value of the mean lift when 



increas ing the wall distan ce; Halll ( 19881) measured consistently positive values for various wall dis- 
tances; IZeng et al. ( 20081 ) obta ined negative valu es (directed towards the wall) in case the sphere 
is placed in the buffer layer. Zeng et al. (2008) note that the classical formulae based on un- 
bounded shear flow fail to predict their DNS results correctly, stating that further investigations 
are requ i red to understand the discrepancy. 

lHalll (|l988l ) showed that the effect of a nearby wall on the lift experienced by a spherical 
body differs significantly depending on the wall being rough or smooth. In particular, it was 
found that the lift significantly decreased when the sphere was positioned in between of spanwise 
oriented, rod-shaped roughness elements. When the sphere was positioned on top of the array of 
wall-mounted rods, however, the measured lift was comparable to the corresponding smooth-wall 
values. 

The difficulties rela ted to the direct measu rement of particle forces as in the studies above 
have been discussed by Muthanna et all (|2005l ). Another, more indirect approach was taken by 



Einstein fc El-Samnil (|l949l ). They approximated the force exerted on hemispheres in an open 
channel flow by pressure measurements on top and near the bottom of the hemispheres. They 
reported positive lift on the hemispheres, and were among the first who stated the relevance of the 
forces on particles in a rough wall to the understanding of sediment erosion. More recent studies 
following this approach present approximatio ns of lift and drag on cubes, spheres and naturally 
shaped stones by local pressure measurem ents ([Holland. Battjes fc Booii . 2005; Hofland fc Battiesl . 
2006t iDetert. Weitbrecht fc JirkaL 1201061 ). These studies have focused on the higher Reynolds 



number regime with particle Reynolds numbers of the order of thousands. 

The investigations discussed so far have for the most part concentrated on the flow around 
single spherical objects. Beyond these studies, a large body of literature exists which deals with 
the characteristics of flow over rough surfaces. Although the precise nature of the fluctuating 
forces acting on individual roughness elements is often not of interest in the context of studies 
of roughness effects, findings from that r esearc h area are of relevance here. A reference for the 
earlier work on roughness is 
numerical studies is given by 



Schlichting] (fl96 5i); a more recent review on the subject, including 
J nncnez "(120041). Some aspects of rough- wall flows at high Reynolds 



number have recently been reviewed by Marusic et al. ( 2010l ). in particular the question whether 
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roughness does indeed modify the turbulence structure in the outer flow or simply provides a 
modified friction velocity. 

The research on rough wall turbulence focuses almost exclusively on the effect of the rough wall 
on the fluid. Some of the key questions of interest are how roughness influences the turbulence 
structure, what are the consequences for scaling, and how can the effect on the fluid be estimated 
from the roughness geometry. 

Numerical studies of rough wall flows are very demanding in terms of computational cost, much 
more so than comparable simulations of flow over smooth walls. Therefore, far less direct numerical 
studies of flow over rough walls have been carried out so far in contrast to flows over smooth wall. 
Recently a direct numerical study of a boundar y layer over surface s rou ghened by rectangular 
spanw ise bars and cubes has been carried out by iLee fc SunsJ (|2007l) and iLee. Sung fc Krogstad 
(|201ir ) . Direct numerical simulation of channel flow over a wa ll similarly ro ughened by spanwise- 



oriented square bar s has b een carried out by iLeonardi et oil J2003, 2007) among others, while 



Orlandi fc Leonardil (|2008) have simulated plane channel flow including different layouts of wall- 



mounted cubes. Direct numerical simulations of ch annel flow with wall vel ocity disturbances 
(actin g as artificial roughness) have been carri ed out by Orlandi et all (120031 ) and Flores fc Jimene2 
(200§). More in line with the present setup, Singh. Sandham fc Williams! ( 2007h have performed 
simulations of open channel flow over spheres in hexagonal arrangement, albeit at considerably 
coarser resolution than the one employed in the present study. 

As a first step to understand the mechanism leading to sediment erosion here we present high- 
fidelity data on the flow over a rough wall with a regular array of fixed spheres. In contrast 
to most previous studies on roughness, the emphasis of the present work is on the effect of the 
turbulence on the spherical elements which form the rough wall, including the characteristics of 
hydrodynamic force and torque. 

The article is structured as follows. In §2] the setup of the simulation is described and basic 
definitions are given. The section also includes a brief discussion of the numerical scheme used. In 
§12 the results are discussed and compared with previous findings in the literature when possible. 
The results of the time and spatially averaged flow field statistics are discussed in £13.11 followed 
by a discussion of the time-averaged three-dimensional flow field statistics in §3.21 The statistics 
of the particle force and particle torque respectively are presented in ^3.31 and ^3.41 jointly with 
their probability density function (pdf) and the local surface distribution of the contribution to 
the mean values. In the discussion the results are related to some degree to forces and torque on 
an area element in a smooth wall channel flow. Conclusions and an outlook are given in §3 



2 Flow configuration 

The flow configuration consists of turbulent open channel flow over a geometrically rough wall. 
The wall is formed by one layer of fixed spheres which are packed in a square arrangement (see 
figure [T]). The distance between the particle centres is D + 2 Ax, where D is the particle diameter 
and Ax is the grid spacing. At y = a rigid wall is located below the layer of spheres. As can be 
seen in figure [T] this rigid wall is roughened by spherical caps that can be defined as the part above 
y = of spheres located at y = D/2 — \/2(D/2 + Ax), staggered in the streamwise and spanwise 
direction with respect to the layer of spheres above. 

The physical and numerical parameters of the simulations are summarised in table [TJ The 
computational domain dimensions are L x / H x L y /H x L z /H = 12 x 1 x 3, in streamwise, wall- 
normal and spanwise direction, respectively. An equidistant Cartesian grid with 3072 x 256 x 768 
grid points is employed. 

One important parameter is the ratio between the domain height, H, and the spheres diameter, 
D. Ideally, a large H/D is desira ble to ensure th at the spheres can be considered as roughness 



and not as obstacles in a channel (jjimenea . 120041 ) . However, from a practical point of view it is 



difficult to reach large values of H/D without increasing excessively the computational cost. In 
this work, two cases are considered: case F10 with H/D = 18.3 and a total of 9216 particles, and 
case F50 with H/D = 5.6 and a total of 1024 particles above the bottom wall. 
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Case U bh /u T Re b Re T D+ D/Ax Ax+ N p T c U bH /H 
F10 15.2 2870 188 10.7 14 0.77 9216 120 
F50 12.2 2880 235 49.3 46 1.07 1024 120 

Table 1: Setup parameters of simulations; U b H is the bulk velocity based on the domain height 
H, U b h is the bulk velocity based on the effective flow depth h defined as h = H — 0.8D, u T is 
the friction velocity, Re b = UhaH/v is the bulk Reynolds number, Re T — u T h/v is the friction 
Reynolds number, D + = Du T /v is the particle diameter in viscous units, D/Ax is the resolution 
of a particle, Ax + is the grid spacing in viscous units, N p is the total number of particles in a 
layer, T c is the time over which statistics were collected. 




Figure 1: Close-up of a section of the computational domain with the geometry of the bottom 
wall consisting of a layer of fixed spheres arranged on a square lattice; (a) case F10; (b) case F50. 



Periodic boundary conditions are applied in streamwise and spanwise directions. At the upper 
boundary a free-slip condition is employed. At the bottom boundary a no-slip condition is applied. 
The spheres are resolved using the immersed boundary method which is described in £ )2.1I 

In order to scale the results, two quantities need to be specified: the friction velocity u T and the 
location of the virtual wall, yo, since for a geometrically rough wall the position of the wall cannot 
be unambiguously defined (cf. iTownsendl . 1971 : Raupach et al . 1991 ). As discussed in detail in 
appendix SjA] we choose to define the position of the virtual wall as yo = Q.8D throughout this 
study. The value of u T is defined by extrapolating the total shear stress T to t = pvd{u) jdy — p(u'v') 
from above the roughness layer (where it varies linearly) down to the location of the virtual wall yo. 
The effective flow depth, h, can be defined as the distance from the virtual wall to the top boundary, 
h — H — yo- The bulk velocity based on the domain height, H, is defined as U b n = 1/H J Q H (u)dy; 

the bulk velocity based on the effective flow depth is defined as U b h = l/h f H (u)dy « U b HH/h. 
Angular brackets are used for the notation of the averaging operator jointly with sub-indexes t,Xi,p 
in order to specify averaging in time, along the direction Xi or over the periodically repeating cells 
of the geometry, respectively. Angular brackets without additional indices refer to quantities 
which are averaged over time as well as spatially over wall-parallel planes, i.e. along the x and 
z directions. The bulk Reynolds number, Re b = l/baH/v, was kept constant at a value of 2870 
and 2880 in cases F10 and F50, respectively. This corresponds to a friction Reynolds number, 
Re T = u T h/v ~ 180 in case of a smooth wall. In the present simulation the value for Re T increases 
to 188 in case F10 and to 235 in case F50. 

The grid resolution of the simulation was in both cases approximately equal to the viscous 
length 5 V — v ju T in all spatial directions. The resolution can therefore be qualified as exceptionally 
fine away from the wall and as reasonably fine in the vicinity of the wall. In the following, 
normalisation with wall units will be denoted by a superscript +. 

The initial turbulent flow field of each simulation was taken from a similar simulation on a 
coarser grid. Subsequently, the simulation was run until the flow reached a statistically-stationary 
state. The simulation was then continued for 120 H/U b n during which flow field statistics as 
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well as particle data such as forces and torques were collected. Entire flow fields jointly with the 
particle data were stored at intervals of about H/UbH- Based on these data a statistical analysis 
has been carried out. If not explicitly stated otherwise, the statistics shown in the following stem 
from the data collected during the runtime of the simulation and are averaged over the entire 
domain including the region within the particles. Details on the different averaging procedures 
used and how they compare are provided in <|B] and SjCj 



2.1 Numerical scheme 



In order to discretise the complex shape of the wall (cf. figure Q]), the pres ent simulations were 



carried out with the aid of a y ariant of the immersed boundary technique (jPeskinl . I1972L [2002) 
proposed by mannl (|2005al ). This method employs a direct forcing approach, where a lo- 
calised volume force term is added to the momentum equations. The additional forcing term is 
explicitly computed at each time step as a function of the no-slip condition at the fixed particle 
surface, without resorting to a feed-back procedure. The necessary interpolation of variable values 
from Eulerian grid positions to particle-related Lagrangian positions (and the inverse operation 
of spreading the computed force te rms back to the Eulerian grid) are performed by means of the 
regularised delta function given bv lRoma. Peskin fc Bergerl (1993). 

A Cartesian grid with uniform isotropic mesh width Ax = Ay = Az is employed which ensures 
that the regularised delta function verifies important identities (such as the conservation of the 
total force and torque during interpolation and spreading). For reasons of efficiency, forcing is 
only applied to the surface of the spheres, leaving the flow field inside the particles to develop 
freely. 

The immersed boundary technique is implemented in a standard fractional-step method for 
the incompressible Navier-Stokes equations. The temporal discretisation is semi-implicit, based 
on the Crank-Nicholson scheme f or the viscous terms and a low-storage three-step Runge-Kutta 
procedure for the non-linear part ( Verzicco fc Orlandi . 19961 ). The spatial operators are evaluated 
by central finite-differences on a staggered grid. The temporal and spatial accuracy of this scheme 
are of second order. 

An important benefit for the present simulation is that the hydrodynamic forces acting upon 
a particle are readily obtained by summing the additional volume forcing term over all discrete 
forcing points. The analogue procedure is applied for the computation of the hydrodyna mic torque . 
T he pres en t nume rical method has been submitted to e xhaustive validatio n tests (jUhlmann 



2003, 2005a. b, 2006a), as well as grid convergence studies ( Uhlmann . 200661 ). In addition, the 



computat ional code has be en applied to the case of vertical plane channel flow with many moving 
particles ((Uhlmann ._ 2008 ) . I n par ticular, this reference also includes a validation against the 



benchmark case of Kim et al. 



(|1987l ). Recently, the present immersed boundary method has been 
successfully implemented and employed in different numerica l codes by other researchers (e.g. 
Lucci. Ferrante fc Elghobashi , 20101 : Lee fc Balachandar , 2010l) . 



3 Results and discussion 
3.1 Flow field statistics 

Figure [5] shows the profiles of the time and plane averaged streamwise velocity component, (it), 
as a function of the vertical coordinate. The results of case F10 and case F50 are compared with 
the reference case S180 of a smooth-wall open-channel flow at i?e& = 2880 and Re T = 183, which 
has been recomputed for the present study. The profiles show the expected effect of roughness 
that is described in various textbooks (|Schlichtingj . (l965t iPopel |2000( ) . As the particle diameter D 



increases, while keeping the value of the bulk Reynolds number Re^ constant, the friction velocity 
increases; the streamwise velocity profile normalised by outer scales flattens (figure [2}i) , and the 
streamwise velocity profile normalised by viscous scales increasingly shifts towards lower values of 
(u) + (figured). Figure [2] shows that a logarithmic layer exists, if at all, only over a small range 
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due to the low Reynolds number considered. The logarithmic law for the flow over a rough wall 
can be written as in the case of a smooth wall with an additional offset AU + that accounts for 
the roughness effect 



(«)" 



1 



In 



y-yo 



A-AU+ 



(1) 



where k and A are constants obtained empirical l y to be n w 0.4 and A s» 5.1 (according to 
experimental findings summarized e.g. in I Jimenea . I2004T) . From the profiles in figure O it appears 
that in case F10 the roughness effect is weak, while in case F50 a stronger roughness effect can be 
seen. 

It is customary to quantify the roughness effect by using the equivalent sand grain roughness 
k s (jSchlichtinsJ . I19361 ). It can be obtained by a fit to the mean velocity profile in the logarithmic 
layer using the following equation 



= C log 10 



y-yo 



B. 



(2) 



wher e B ps 8.48 and C = 1/k In (10) ps 5.75 are empirically obtained values fcf. lShockling. Allen fe Smits 
2006). At high enough Reynolds numbers k s becomes constant, i.e. lim^e^oa k s = k soo , defining 
the so-called fully rough flow regime. The specific value of k soo is a property of the surface that de- 
pends on the characteristics of the roughness, such as shape, arrangement or roughness area ratio. 
Flow over roughness can be classified as hydraulical ly smooth, transitionally rough or fully rough 
according to a small, moderate or high value of . iNikuradsd (|l933[ ) gave values of 5 < kf^ < 70 
to define the transitionally rough flow regime. However, these values should be taken with care as 
th e transition migh t be influ e nced by the specific characte ristics of the roughness (cf. discussion 



Bradshaw , 120001: J lmene 3. l2004t ISnockling et all 2006, among others). In particular, it has 



been speculated that a uniformly sized, structured arrangement of roughness elements as in the 
present case might lead to a sharp transition from hydraulically smooth to the fully rough regime 
(|Colebrookl Il939t iJimenea . I2004T ). Figure [3] shows the transition from the hydraulically smooth 
flow regime to the fully rough flow regime obtained in different experiments. It shows the offset 
AU + as a function of k^. At low values of k^ the effect of roughness should be negligible and 
correspondingly AU + approaches zero. In the fully rough flow regime the roughness effect should 
purely depend on k^. Comparing equations JT]) and ([2|) a formula for AU + can be derived for 
the fully rough regime, 



AU+ = C log 10 (fc+J-B + A, 



(3) 



with the constants A, B and C as above. The relation ([3]) above is shown in figure [3] jointly with 
results from experiments. 

For the present simulations the values of AU + can be obtained by the vertical shift of the 
mean velocity profiles in the log- region of figure &&b)- They are 1.03 and 4.85 for cases F10 
and F50, res pectively. Howeve r, the value of k so? /D fo r the present arrangement of s pheres 
is unknown. ISchlichtind (|l936l ). iLigrani fc Moffatl (|l986l ) and IPimenta. Moffat fc Kavsl (|l975l ) 
found a value of k soo /D ~ 0.6 3 for flow o v er sp heres in hexagonal pac king, while somewhat 
larger values were obtained by ISingh et all 1)20071 ) (k soa /D = 0.77) and iDetert et all (|2010al ) 
(k s /D — 0.81). For flow over spheres in r andom packing the valu e s gen e rally obtain e d var y 
in the range of k soo /D = 0.55 to 0.85 (cf. iMufioz Goma fc Gelharl. Il968t iGrass et all Il99ll) . 
Few studies exist that use values of k soo /D = 1 for structured ( Einstein fc El-Samnil . Il949h or 
random arrangements ( Nakagawa fc Nezul . Il977t ). Figure [3] shows the pair of values (kf^, AU + ) 
for cases F10 and F50 when approximating k soo / D by the value found for a hexagonal packing, i.e. 
ksoo/D = 0.63. The error-bars indicate the range of k soo /D = 0.55 to 1 as found in the literature. 
It can be seen that case F10 approaches the hydraulically smooth flow regime while case F50 is in 
the transitionally rough flow regime. 

A similar conclusion can be reached by analysing the profiles of the root-mean-square of the ve- 
locity fluctuations normalised with u T that are shown in figure[4ja). In case F10, the profiles of the 
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(u)/u bh {y-yo)/s v 

Figure 2: Time and spatially averaged streamwise velocity component (u) of case F10 (dashed line) 
and case F50 (dashed dotted line) in comparison with smooth wall open channel flow (continuous 
line); (a): normalised with U^h as a function of (y—ya)/h; (b): in semi-logarithmic scale normalised 
by 8 V and u T ; the position of the particles top are marked with horizontal (a) and vertical (6) solid 
lines. 




Figure 3: Roughness function for se veral transitiona lly rough surfaces as a function of the Reynolds 
num ber based on kt^ adap ted from l Jimenez! ( 2004) . o. lNikuradse ( 19331) . uniform sand, pipe flow ; 



V. iLigrani fc Moffatl (|1986), uniform densely-packed spheres, boundary layer; IShockling et al. 
(|2006f ). honed aluminium, pipe flow; ■, present simulation with k soo /D = 0.63, error -bars show 
the ra nge of k soo /D = 0.55 to 1; dashed line, AU + — 5.75 log (1 + 0.26/c^) proposed bv lColebrook 
(|l939t ); dashed dotted line, relation © with A = 5.1, B = 8.48 and C = 5.75. 
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three velocity components almost collapse with the smooth-wall results, indicating that, indeed, 
the flow over the relatively small roughness elements in this case can be considered as nearly hy- 
draulically smooth. In case F50, some differences with respect to the smooth wall case are evident. 
The near-wall peak in the streamwise fluctuation profile decreases but it is still visible. This indi- 
cates that the flow is in the transitionally rough flow regime si nce experiment s in the fully rough 
flow regime present no clear peak (see for example figure 5 of I Jimenez! 120041 ) . The wall- normal 
and spanwise fluctuations present slightly higher values near the wall than the corresponding ones 
in the smooth- wall case. Therefore, the anisotropy of the fluctuations near the wall is smaller than 
in the smooth- wall case. This tendency of roughness to make the fluctuations more isotropic is 
a phenomenon which has been o ften reported in the literature (e.g. IPoggv. Porporato &: Ridolfi" 



2003; Orlandi fc Leonardi . 20081 ). Also in case F50, above (y — yo)/h ~ 0.4 all three components 



ag ree well with the values of t he smooth-wall case 



Orlandi fc Leonardil ([20081 ) discuss the velocity shift AC/+ as a function of v rms at the roughness 
crest. The present simulations result in pairs (v rms , AU + ) of (0.10, 1.03) (0.46, 4.85) for case 
F10 and F50 respectively which agree well within the scatter of the reported data (graph not 
shown) . 

Figure \M,b) shows the profiles of the Reynolds stress, (u'v'), normalised by v? T . The Reynolds 
stress profile (u'v') of case F10 nearly collapses with the profile of the smooth-wall case. In case 
F50 a slight increase and a small shift towards the wall of the near-wall peak can be seen which 
could be an effect of the higher value of Re T in this case. 

In order to study the near wall behaviour of the velocity fluctuations, a close-up of the profiles 
shown in figure Ufa), is plotted as a function of (y — yo)/D in figure [S] Additionally, the profiles of 
the root-mean-square of the pressure fluctuations, p rms / (pv%), are included. Note that in contrast 
to figure S] the profiles shown in figure [S] are obtained from snapshots of the flow field and obtained 
by averaging over cells outside of the particles as described in detail in S|B]and $Cj The amplitudes 
of the fluctuations of the three velocity components present similar values below the virtual wall, 
(y — yo)/D < 0. These are much smaller than the values above the virtual wall, and they are 
somewhat larger in case F50 (u l rms ~ 0.1u T ) compared to F10 (u z rms < 0.05u r ). On the contrary, 
the pressure fluctuations within the roughness layer for both cases present values which are similar 
to the values above the roughness layer. 

R ecall that also in the case of a smoot h wall the pressure fluctuations are non-zero at the 



wall dKim. Moin fc Moserl Il987t iKiml . Il989h . Near the top of the roughness elements, i.e. around 



(y — yo)/D = 0.2, the profiles of p rm s exhibit a peak which is barely visible in case F10 and more 
pronounced in case F50. In case F50 the value of the peak is higher by a factor of two compared to 
case F10. Note that the peak is, to some extent, a consequence of the three-dimensionality of the 
time-averaged flow field around the particle; this point is further elaborated in tyd.'2{ The pressure 
fluctuation profiles of case F10 and F50 (when plotted as a function of (y — yo)/h) approach each 
other with increasing wall distance (not shown) . They converge to the profile obtained in the case 
of a smooth wall in the outer part of the flow. 

3.2 Three-dimensional time-averaged flow field distribution 

Since the geometry of the roughness is three-dimensional the time-averaged flow field in the near 
wall region also varies in all three directions. In the following some characteristics of the time- 
averaged flow field obtained from 90 snapshots are discussed. In addition to the averaging in time 
the fields were averaged over periodically repeating boxes centred on the particles (henceforth 
indicated by the symbol (-)tp)- For simplicity this will be simply referred to as time averaging 
below. 

Figure |5] shows the distribution of the three-dimensional time-averaged streamwise velocity, 
(u)t P , for both cases. Two different (x, y)-planes are shown. The first one contains the centre of 
the particles of one streamwise row (figure |6ji and figure [6^). The second one is located in between 
two streamwise rows of particles (figure [Sf) and figure [SJi). As can be expected the flow field is 
very different from a single sphere in an unboun ded turbulent flow or in a channel close to a wall 
(jBagchi fc Balachandarl 120041: IZeng et all 120081) . 
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Figure 4: (a) Root-mean-square of velocity fluctuations of case F10 and case F50 normalised by 
u T in comparison with results of smooth wall as a function of wall distance; curves from left to 
right wall-normal (v rms /u T ), spanwise (w rms /u T ) and streamwise (u rms /u T ); (b): distribution of 
Reynolds shear stress (u'v') as a function of wall distance. Legend as in figure [2] The position of 
the particle tops are marked with horizontal solid lines. In (b) the straight dashed line is included 
to guide the eye. 




Figure 5: Root-mean-square of velocity and pressure fluctuations of case F10 and case F50, nor- 
malised by u T and pu^ respectively, as a function of (y — yo)/D; averaging has been carried out 
over fluid cells only, for details see $B] □ : u rms /u T , V : v rms /u T , A : w rms /u T , and O : Prms/ pu 2 ] 
solid lines and empty symbols: case F10, dashed lines and full symbols: case F50; horizontal line: 
position of particle tops in both cases. 
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The sheltering effect of the neighbouring particles causes the flow velocity to decrease rapidly 
close to the roughness tops and leads to marginal flow velocities within the roughness layer. The 
highest velocity gradients are prod uced in the vicinity of the roughness tops. Similar observations 
were made in the experiments of IPokraiac &: Manes (2009) who studied a comparable particle 



arrangement at bulk Reynolds numbers of order 10 and a ratio of h/D ~ 3.5. Also similar to 
their results is the formation of a recirculation between two spanwise rows of spheres that extends 
over the entire spanwise direction. The shape of the recirculation is similar in both of our present 
cases, however the strength differs. In case F50 the backflow velocities reach values as low as 
(u)tp ~ — 0.4u T in figurc^c) and (u) tp ~ —0.2u T in figure^d). In case F10 the magnitude of the 
backflow velocity is below 0.05w T . The recirculation can also be observed in figure [7] that shows 
streamlines of the mean flow projected into the same planes shown in figure[Sl i.e. by computing the 
streamlines using only (u} tp and (v)t P , together with contours of the time- averaged pressure field. 
The pressure distributions in both simulations are similar. However, in case F10 the magnitude 
of the pressure, (p) t p / (pv%) , is a factor of two smaller compared to case F50. Please note that 
the three-dimensional time-averaged flow is not fully converged and at the location of the planes 
shown in figure [7j there is a weak net flow in the spanwise direction with a maximum amplitude 
of (w) tp w 4 ■ 10~ 4 [/ fc/l (6 • 10~ 4 [/ 6/l ) in case F10 (F50). This net flow is within the range of the 
statistical uncertainty. 

A question of interest is how far the three-dimensionality of the bottom wall directly influ- 
ences the flow. Figure [7J already shows that one particle diameter above the roughness tops 
the time-averaged pressure field is still visibly affected. In order to quantify the effect of three- 
dimensionality, the difference between the time-average of a field (4>)t P (where 4> can stand for 
either pressure or one of the velocity components) and its time and plane-averaged value, (</>), can 
be defined, viz. 

4>" = {4>) tp - (0) . (4) 

Note that the quantity defined by ( PQ) is sometimes cal led spatial disturbance in the context of 



the double-averaging methodology (jNikora et all . 120011 1. From equation (j3J the corresponding 



standard deviation 0"ms = V { ( t ) " ( t ) ") can be computed. In both cases, F10 and F50, the standard 
deviation of pressure p" ms drops by several orders of magnitude in between y — D and y = 2D, 
as shown in figure |SJ The same is true for the velocity field (not shown). Therefore, the time- 
averaged flow statistics appear to be essentially one dimensional beyond wall distances of 2D. 
This is s omewhat smalle r than the values reported in previous investigations of flow over rough 
walls (cf. I Jimenez , 2004 ) which might be related to the low values of D + and Re T considered. 



3.3 Statistics of particle forces 

The hydrodynamic force, F, acting on a particle is defined as 

r-ndr- f p tot n dT, (5) 
r Jr 

where T is the sphere's surface, n, is the surface normal vector, r = pv {djUi + diUj) is the 
viscous stress tensor and p tot is the pressure. The latter can be split into two parts p tot = p + pi, 
where pi represents the linear variation in streamwise direction which results from the imposed 
pressure-gradient that drives the flow, and p corresponds to the three-dimensional instantaneous 
fluctuation. The first term on the right hand side of equation ((5|) is the force due to viscous 
stresses, the second term is the force due to pressure. A sketch that illustrates the definition of 
the force on a particle can be seen in figure IUa) . 

In order to scale the hydrodynamic forces, reference quantities need to be defined. For the 
present case of particles within a roughness layer the s ubject is a matter of discussion and several 
definitions have been proposed in the literature (see Hofland et all [2005). Here, the reference 



force is defined as Fr — pu^An with the reference area Ar = L X L Z /N P . 

Table [2] summarises the particle force statistics of the two cases, where Cp is the mean force on 
a particle in a^-direction normalised by Fr. As can be seen, the mean values of the forces acting in 
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Figure 6: Distribution of the time-averaged streamwise velocity in a periodic cell, (u} tp , normalised 
by u T in x-y planes; (a, c) show a plane through particle centres, (6, d) show a plane centred between 
spheres. Dashed lines: iso-contours of (u)t P at values of -0.5 to 2.1 in steps of 0.2; black solid line: 
streamwise velocity contour at — 10~ 3 . Panels (a) and (b) show case F10, (c) and (d) show case 
F50. The direction of the bulk velocity is from left to right in all panels. 



Case Cp C V F Cp a F ap/F R <r v F /F R a z F /F R Sp S y F S Z F K X F K V F K Z F 
F10 1.04 0.19 0.00 11° 0.57 0.20 0.66 0.18 1.80 0.01 10.13 19.08 9.92 
F50 1.15 0.37 0.00 18° 1.32 0.66 1.26 0.06 0.26 0.01 4.98 5.68 4.29 

Table 2: Statistics of particle forces in case F10 and case F50, where Cp* — (F Xi /F R ) is the 
normalised mean force component in the Xi-direction, a — arctan(Cj£/Cf.) is the angle of the 
resulting force with respect to the ai-axis, u x F is the normalised standard deviation of the force in 
Xi, Sp 1 and Kp 1 are the skewness and kurtosis of the respective force component. 
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Figure 7: As hgure|6]but showing the time-averaged pressure held and corresponding streamlines; 
dashed lines: iso-contour lines of pressure, (p)t P / (pv^.), from values of -4 to 4 in steps of 0.4; 
continuous lines: streamlines in the plane computed from (u)t P and (v)t P - The direction of the 
bulk velocity is from left to right in all panels. 
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Figure 9: Sketch illustrating the definition of force ( — ►) and torque ( — ►►) on a particle (a) 
and a square surface element in a smooth-wall channel (b) 
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the streamwise direction (henceforth also called "drag" ) and the wall- normal direction ( "lift" ) are 
positive. Since the mean forces are directly related to the mean flow through the time-averaged 
version of equation ([3]), it is possible to shed some light onto the mechanisms that lead to drag 
and lift by analysing figures [5] and [7J A more detailed picture can be obtained from figure [TU] and 
figure [TT1 which show the distribution on the sphere's surface of the stress leading to drag, td, and 
lift, tl, viz. 

r D = ((T) t .n-(O t n)- ei , (6) 
r L = «T) t -n-<p to V)-e 2 , (7) 

where ej is the unit vector in the Xi-direction. The stresses in figures [TU] and [TTJ are normalised 
by Fji/Agph, where A sp h = ttD 2 is the surface area of the sphere; by virtue of this normalisation 
the total integral of the quantities shown in the figures yields the force coefficients Cp given in 
table [3] Please note that the results of case F50 appear less smooth due to the smaller number of 
particles, and therefore a smaller number of samples. 

Figure [UJ shows that the local stress contributing to drag is similarly distributed over the 
particle surface in both of our present flow cases F10 and F50. One can observe a region of strong 
positive values with the largest magnitude centred around a position slightly upstream of the 
particle tops. From figure ITUT a. c) we can see that this region of high positive local contributions 
to drag is slightly elongated in the spanwise direction. It results from the wall-normal gradients 
of the average streamwise velocity component which are particularly important in the upper part 
of the sphere as well as from the high pressure values found near the upstream side of each sphere 
(cf. figures [()] and [7]) . On the downstream side of the particles, still in the upper hemisphere, a 
smaller region with weak negative contributions to drag is found, as a result of the recirculation 
region. In most of the lower (near-wall) half of the spheres, the contour lines of the local drag 
contribution are roughly oriented in the wall-normal direction, changing sign slightly downstream 
of the cross-stream plane passing through the particle centre. In this context it should be noted, 
that the driving pressure gradient dpi /dx < makes a weak but non- negligible contribution to the 
drag which can be quantified as approximately 2% (9%) of Cf. in case F10 (F50). Therefore, non- 
negligible values of local contributions to drag are expected even in relatively quiescent regions, 
as is the case inside the roughness layer. 

The qualitative and quantitative similarity of the distribution of tjj in both cases F10 and F50 
results in similar values for the drag coefficient in both cases (cf. tabled. In particular, the overall 
drag coefficient Cf. in case F10 is close to unity, increasing to 1.15 in case F50. These values are 
a result of the weak contribution of the drag on the rigid wall below the layer of spheres to the 
total drag on the wall and the choice of the reference force. The drag coefficient as defined in the 
present study can be approximated as 

o- ~ v l tot) + v ^ (8 ) 
° F /^T~ ' (8) 

where vj tot ^ is the total volume occupied by fluid in a periodic cell around a particle, V sp h is 
the volume occupied by a particle. The approximation (|8]) neglects the streamwise component of 
the shear force acting on the bottom wall in addition to the drag due to the periodic part of the 
pressure acting on the spherical caps. Evaluating this geometrical relation ([8j yields 1.04 (1.15) 
for case F10 (F50). 

Positive values for the lift coefficient, as observed in the present simulations (cf. tabled]), can 
be explained by two mechanisms. The approaching flow accelerates in the frontal part until the 
top of the sphere and from then on it decelerates. This fact is reflected in the curvature of the 
streamlines (figures [UJi, c), yielding a pressure distribution which exhibits lower values of pressure 
near the particles tops (figures 02, c), and therefore a positive lift. In addition to pressure, shear 
might lead to a positive lift. As can be seen in figure [UJ the flow field above the spheres is 
asymmetric (with respect to a cross-sectional plane through the particle centres) as a result of 
the recirculation behind the particles. Therefore, the friction on the upstream side of the particle 
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(in the upper hemisphere) is expected to be higher compared to the corresponding friction on the 
downstream side, contributing positively to the lift on the particle. The pressure differences as 
well as the asymmetry of the flow seem to be more pronounced in case F50 than in case F10 and 
might explain the observed increase in the lift coefficient. 

Figure [TT] shows the distribution on the sphere's surface of the stress leading to lift, tl. The 
shape of the contours is again similar in both cases, however, the magnitude of the stress tl seems 
to be significantly larger in case F50, leading after integration to the factor of two presented in 
table [21 The spatial distribution is characterised by one dominant patch of each positive and 
negative values of T£, the maximum of both being located on the (x,y) symmetry plane, the 
former (positive) near the particle top, the latter (negative) shifted upstream by approximately 
45° (30°) in case F10 (F50). From the contours, it appears that the flow below the virtual wall 
contributes little to the lift. 

In order to quantify the contribution integrated from the bottom of a particle up to a certain 
fraction of its diameter, we can define a cumulative function 

s*(y) = ^ f dfldy, (9) 

z Jo Jo 

where T^{y,ff) stands for either T£ or tjj evaluated at a position on the sphere's surface given by 
the wall-distance y and an azimuthal angle 9 in the wall-parallel plane. Figure [T^] shows Sl and 
Sd normalised by the net values of lift and drag, respectively. The contribution to the net drag by 
the flow in the lower half of the sphere is small in both cases, the cumulative drag value increasing 
monotonically and with increasing slope from the wall to the top of the sphere. Conversely, the 
cumulative contribution to the lift first increases with increasing wall-distance up to values of 
approximately 25% (40%) of the total in case F10 (F50) at y ~ 0.5/?, before decreasing again to a 
small value at y « 0.9D. Beyond that, in a small area surrounding the top of the sphere, is where 
most of the net lift is generated. In case F50, the lift increases with respect to case F10 more than 
the drag, which leads to a higher angle, a = arctan(Cj^/Cf.), of the resulting force (cf. table [2]). 

The spanwise force should be zero for symmetry reasons. In both cases the calculated mean 
spanwise force coefficient is more than two orders of magnitude lower than the drag coefficient. 
This fact provides confidence in the convergence of the statistics. 

Additional support to the mean forces just discussed is provide d by compar ison to experimental 
measurements performed in a somewhat similar configuration by Hall (1988). In that study, the 
mean lift on a particle near a boundary was measured in a wind tunnel with smooth as well as 
rough walls. The interesting case for the present discussion consisted of a sphere of diameter D 
placed in-between spanwise rods of diameter D r evenly spaced out with a distance D r . Figure 
[TBI presents the comparison of the mea n lift norma lised by pv 2 as a function of D + , between the 
values obtained in the experiments of Hall (|l988l) and the present simulations. In spite of the 
different setups, the lift obtained in case F10 is perfectly consistent with the measurements while 
the lift obtained in case F50 is somewhat lower. The reason for this might be that in the setup of 
the simulations the neighbouring spheres are closer producing an increased sheltering effect. This 
is also supported by the ex perimental observation that lower lift values are obtained when the 
value of D r /D is increased (jHalll . 1 19881 ). 

In contrast to the direct relation between the mean flow field and the mean forces on a particle, 
a similar straightforward relation between the statistics of the fluid velocity fluctuations and of the 
particle force fluctuations cannot be derived from equation ([5]). This is due to two factors. First, 
the definition of the standard deviation (and higher order moments) of the force fluctuations is 
non-linear. Second, the integrals in equation §5§ act like a filter in the sense that not all scales 
participate in creating force fluctuations on a particle. For example, flow scales much smaller than 
D might cancel out in the integral sense as will be discussed in detail in §3.41 In spite of this 
observation, a direct relation between flow velocity statistics above the bed or behind an obstacle 
is often ass umed in the literature in order to estim ate the intensity of force fluctuations on a 
particle (cf. [Papanicolaou et al , 2002 ; Garcial [20081 ). 

In the present simulations we observe that the standard deviations for the streamwise and 
spanwise components of the particle forces are of similar magnitude in both cases F10 and F50 
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Figure 10: Spatial distribution of td, normalised by Fn/A sp h- x and z are the coordinates with 
respect to the particle centre. The contour lines shown correspond to [-0.5, 0.0, 0.5, 3.5, 7.5]; 
Panels (a) and (b) show case F10, panels (c) and (d) show case F50. 
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Figure 11: As figure [TU1 but for tl- The contours lines shown in (a) and (b) are at values from 
-2.5 to 2.5 in steps of 1, in (c) and (d) at values from -7.5 to 10.5 in steps of 3; in (d) additionally 
the contours at the values -0.5 and 0.5 are shown. 
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Figure 12: Cumulative function of the stress contribution to the mean value of drag, lift and 
spanwise torque on a particle as a function of y and normalised by its maximum value. — □— , 
drag; — • — V — ' — > lift; • • • A • • • , spanwise torque. Panel (a) shows case F10, panel (b) shows 
case F50. 
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Figure 13: Comparison of (F y ) / (pv 2 ) as a function of D + of case F10 and case F50 (sol id symbols) 
with mean lift on a sphere placed in between roughness elements in a boundary layer bv lHaill (1988) 
(open symbols); present simulations: ■; experiments lHalll (|l988 ): V : D r — 5/31?, : D r = D, 
<: D r = 2/31?, where D r is the radius of the rods spaced with D r upstream and downstream of 
the sphere. 
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Figure 14: Normalised probability density functions of force fluctuations, (a): case F10; (b): case 
F50; continuous line: Gaussian distribution, dashed line: F^/ap; dash-dotted line: F' y ju F \ dotted 
line: F'Jap. 



(cf. table [2]) . It is also found that the standard deviation of lift in both cases is roughly half 
the value of the other two components. Overall the intensity of the fluctuations in case F50 is 
more than a factor of two larger than in case F10. Thus, the particle force fluctuations in the 
present case do not seem to scale directly with the intensity of the plane- and time-averaged fluid 
velocity fluctuations (cf. figure |4|), since u rms is larger than w rms over most of the flow depth, and 
especially close to the wall. Furthermore, the difference in the fluid velocity fluctuation intensities 
between case F10 and F50 is very small compared to the above stated difference in the particle 
force fluctuation intensities. It can therefore be concluded that a direct link between fluid and 
particle force fluctuation intensities cannot be inferred in the present cases. 

The results of skewness and kurtosis of the force distributions (cf. table [2} are now discussed 
jointly with the probability density function (pdf) of the particle force fluctuations shown in figure 
IT41 For both cases F10 and F50 the highest skewness is obtained for the lift, i.e. Sp. In other 
words, large positive lift fluctuations are significantly more likely to occur than large negative lift 
fluctuations. This is clearly visible in figure H4Ta). where lift events of several standard deviations 
higher than the mean have a non-negligible probability of occurrence. In case F50, this effect is 
not as strong as in case F10 (cf. figure [T4b). and accordingly the value of the skewness S V F is lower 
in the former case. The small positive skewness of the drag indicates similarly that instantaneous 
high drag events are more likely compared to low drag events. For this component, however, the 
effect appears to be much weaker as compared to lift. Finally, symmetry arguments again lead to 
the conclusion that Sp should be zero, and this is indeed the case. 

The kurtosis of all profiles is rather large indicating a strong intermittency of the forces, i.e. the 
pdfs in figure [14] exhibit much longer tails than a Gaussian distribution. However, as the spheres 
become larger the values of skewness and kurtosis approach the Gaussian values of zero and three. 
This trend might be due to the fact that the force on the particle is an integral quantity, and 
as mentioned before, small intermittent events might be averaged out. This argument is further 
elaborated in §3.41 below. 

Th e present results might be compared to the experimental data provided by lMollinger fc Nieuwstadt 
(|1996T ) for lift fluctuations on a single sphere with D + = 2.9 positioned on top of a smooth wall. 
Although their flow configuration is somewhat different (no sheltering effect, turbulent boundary 
layer) they also report positive values for the skewness (Sp — 1.2) and high values of flatness 
(Kp = 7.0). Furthermore, the pdf of the lift fluctuations in their study is of similar shape to the 
one obtained in the present case F10. This qualitative agreement suggests that the present results 
might be relevant to a broader range of flow configurations, e.g. different sphere arrangements or 
packing densities. 
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Case C T C V T Cf, o T /T R (J y T /T R <t t /T r S t S v t S t K t K v t K t 
F10 0.00 0.00 -0.98 0.21 0.04 0.36 0.01 -0.01 -1.04 6.46 6.17 4.72 
F50 0.00 0.00 -0.73 0.17 0.11 0.27 -0.01 -0.01 -0.76 3.75 4.91 3.37 

Table 3: Statistical moments of torque on particles in case F10 and F50. = (T Xi )/T R is the 
normalised mean torque component in the a^-direction, a T * is the standard deviation of the torque 
in Xi direction, S T * and K T * are the skewness and kurtosis of the respective torque component. 



3.4 Statistics of particle torque 

The hydrodynamic torque T acting on a spherical particle with respect to its centre is defined as 
follows: 

T = J r c x (r ■ n) dT , (10) 

where r c = [x c , y c , z c ) is the distance vector from the particle centre to an element of the surface 
r. It should be noted that - contrary to the definition of the total particle force © - the pressure 
does not enter the integral (fTTJ)) . since in the present case the differential pressure force —p tot nds 
acting on a surface element ds, is always directed towards the particle centre. Based on the 
reference force F R given in ^3.31 the reference torque is defined as T R = F R r Rl where r R is the 
distance from the particle centre to the virtual wall, r R = Uq — D/2. The quantity T R will be used 
in the following for the normalisation of the various torque-related statistical values. A sketch 
that illustrates the definition of the torque on a particle can be seen in figure [5fa). 

Table [3] shows the statistical moments of the torque acting on the particles. Here C T l is the 
mean torque in the a^-direction normalised by T R . Once more, due to symmetry the only non-zero 
component of the mean torque is expected to be C T . The table shows that negative mean values 
for the spanwise component are obtained. These negative values of C T are expected for the torque 
on a particle in positive shear (cf. figure [5] and figure ■ The torque coefficient as it is defined 
above takes values close to —1 for case F10, while it is approximately 25% lower in magnitude in 
case F50. 

In order to analyse these integral results in more detail figure [15] shows the distribution on the 
sphere's surface of the stresses leading to spanwise torque, 

tt = t l x c - T D y c . (11) 

The distribution of tt in figure [15] is in both cases similar in shape and values of the contours. 
As for the distribution of Td (cf ■ figure [TU]) , a shift towards the particle front can be observed 
for the minimum values of tt near the particle top. This shift is more pronounced in case F50. 
Negative values of tt occur almost exclusively in the upper part of the particle. Thus over most 
part of the sphere tt is positive, but low in magnitude. The cumulative contribution function of 
tt, denoted by St (cf. equation [5]), which is also shown in figure [T2l reveals that when integrating 
the contribution of tt in the lower part of the sphere it adds up to approximately — 0.15Cf. in 
the vicinity of the virtual wall. In both cases the values of tt are predominantly negative for 
wall-distances above y ~ 0.8D, such that St vanishes around y — Q.9D. It can therefore be 
argued that the net spanwise torque C T is generated in the surface area between a wall distance 
of y = 0.9D and the particle top (i.e. the region highlighted by a dashed line in figure [TB]). 

Before turning to the discussion of the torque fluctuations, first a simple model is introduced, 
which allows us to elucidate some of the characteristics of the torque fluctuations by considering 
the scales of flow motion that lead to the generation of torque on a particle. It should be noted 
that other authors have previously investigated the relation bet ween flow struc tures at different 
scales and the forces/torque exerted upon sediment particles (e.g. Hoflandl . 2005). Here we employ 



a somewhat different approach which allows us to use data from a smooth- wall flow. In particular, 
we first analyse drag and torque fluctuations experienced by a square wall-element in channel flow 
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Figure 15: Spatial distribution of tt, normalised by Tji/A 8p h] x and z are the coordinates with 
respect to the particle centre. Contour lines are shown at values of [-15, -9, -3, 1] in all plots; in 
(b) and (d) additionally the contour line at zero value is shown. The dashed line indicates the 
location of y = 0.9-D. Panels (a) and (b) show case F10, panels (c) and (d) show case F50. 
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Figure 16: Statistical moments (a showing the root-mean-square value, b the kurtosis) of force 
and torque obtained from the simple model described in the text and of actual values of the 
corresponding torque obtained on the particles in case F10 and F50. In the former case (smooth 
wall) the integration is performed over a square wall element of side length s, the open symbols and 

lines corresponding to: □ , T x \ — • — V > Ty\ • • • A« • • , .7>. The forces are normalised 

by pu 2 s 2 , torque is normalised by l/2pu 2 s 3 . In the latter case (rough wall) the integration is 
taken over the particle surface (as given in equation [TO)) with the filled symbols corresponding to: 
▲ , T x ; ▼ , T y ; M. T z . Note that in cases F10 and F50 the reference length s is taken as ^/Ar. 
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with a geometrically smooth wall, systematically varying the linear dimension of the wall-element. 
Subsequently, the obtained statistical results are related to the corresponding statistics of the 
components of the torque acting on a spherical particle in our main simulations. 

Analogously to the definition of force and torque on a particle the force in the x and z direction 
on a square surface element in a smooth-wall channel with area A s = s 2 can be defined as 



r+s/2 r + s/2 r+s/2 r + s/2 

Fx= I I t xv \ Q dxdz, Fz= / T zy \ dxdz, (12) 

J -a/2 J-s/2 J -a/2 J -a/2 

where s is the side length of the element and ry| _ are the components of the stress tensor at 
the wall. The torque on the element with respect to its centre can be defined as 



f r s/ i \ 

T v = / [r s x T zy \ y=0 -r s z T xy \ y=Q J dxdz, (13) 



/2 ,a/2 
-a/2 J-s/2 

where r s is the direction vector with respect to the centre of the area element. A sketch that 
illustrates the definition of the force and torque on a square element in a smooth- wall channel can 
be seen in figure M,b)- 

Drag and spanwise force on the smooth-wall element are expected to be mostly affected by 
velocity scales in streamwise and spanwise direction, respectively, that are of sizes similar or larger 
than s. The effect of velocity fluctuations at length scales much smaller than s will tend to cancel 
out due to the integral character of the force (fT2"]) . Thus the highest value of force fluctuation 
should be expected for smallest values of s, as the contribution of the smaller scales is lost for 
larger values of s. Conversely, due to the cross-product in (p~3|) the torque on a smooth- wall surface 
element is mostly affected by wall normal vortical motions of sizes comparable to s. The effect of 
much smaller and much larger scales will cancel out or lead to only small values of torque. Thus 
for small as well as high values of s small values of T y are expected. At some intermediate value of 
s, the characteristics of wall normal vortical motions should be most efficient in generating torque, 
leading to maximum values of Ty. Figure [THT a) supports that hypothesis. It shows the normalised 
standard deviation of the forces and torque on the surface element, ojt, cjr and Oj- normalised by 
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pu^s 2 and l/2pu 2 s 3 , respectively. On a square element with s + ~ 70 torque appears to be most 
efficiently produced. This value is somewhat larger than the average distance between the low 
speed and high speed streak close to the wall which is commonly found to be of order b0v/u T . As 
can be seen in figure the kurtosis of the above quantities monotonically decreases with the 

size of the surface element indicating that the intermittency of the small scales is larger than that 
of the large scales. 

A direct analogy between this smooth-wall model and the force and torque on a particle is 
not fully justified, as the flow and the geometry are more complex in the present case. However, 
some of the characteristics of the particle torque statistics obtained for the present cases can be 
explained with the aid of such a simple model as will be discussed in the following. Please note that 
only one torque component can be defined for a plane wall element (here T y ), in addition to the 
two in-plane forces considered (T x and J- z ). A correspondence with the three torque components 
acting on a spherical particle is established when considering the plane wall element as being 
located at the top (i.e. the pole located at y — D) of the particles in case F10 and F50. The 
componentwise correspondence is then: J- x — > — T z , T y — > T y , J- z — > — T x (cf. figure [9]). 

The normalised standard deviations of the particle torque components shown in table [3] are all 
non-zero as can be expected. The amplitudes of the fluctuations of the streamwise and spanwise 
torque components are found to be the largest, while the wall- normal component is significantly 
weaker. Compared to the small-sphere case (F10), the streamwise and spanwise components are 
smaller in the large sphere case (F50), by 20% (ct|,/Tr) and 25% (ct^/Tr), respectively. Contrarily, 
the wall- normal value ct^/Tr is significantly larger in case F50 than in case F10 (nearly by a factor 
of three). 

Figure which has already been partially discussed above, also shows the second and fourth 
statistical moments of particle torque fluctuations as a function of particle size. As can be seen the 
standard deviation (figureQljJi) of the wall normal torque acting on the particles, T y , matches rather 
well the standard deviation of the wall normal torque exerted on a comparable-size square element 
in the reference smooth- wall flow, T y - In addition, the figure shows that the standard deviation 
of the spanwise torque, T z , compares very well to the standard deviation of the drag exerted on 
a square element in the smooth wall case, J- x . Concerning the streamwise component of particle 
torque, T x , it is found that its standard deviation is somewhat larger than the standard deviation 
of the spanwise force fluctuations in the smooth wall model, T z \ however, both exhibit a similar 
decreasing trend with increasing values of the length scale. The overall good agreement between 
fluctuation intensity of forces/torque acting on an element of a smooth wall and the corresponding 
torque components of the particle in case F10 and case F50 is interesting for several reasons. 
First, it suggests that the significant torque fluctuations are generated in a rather limited region 
around the particle tops where apparently to some extent the analogy with the hydrodynamic 
action on a wall-parallel square element holds. In particular, the present simulations F10 and F50 
provide two data points in the hydraulically smooth and transitionally rough flow regime, which are 
fully consistent with the existence of a length scale/particle size of maximum wall- normal torque 
generation, as suggested by the simplified model. Secondly, if the above analogy is accepted, then 
it implies that the response of the particle torque fluctuations to the near-wall turbulent flow can 
indeed be described as a selective filtering effect, mainly characterised by a single length scale (the 
particle diameter). 

Normalised pdfs of the particle torque fluctuations are shown in figure [17] It can be seen 
that the curves for all three torque components in both cases F10 and F50 approximately match 
the curves of the corresponding force/torque components of the smooth- wall model (evaluated 
with a side- length s matching the respective length V ' Ar), thereby further corroborating the 
analogy. Concerning the shape of the particle torque pdfs themselves, it is observed that the two 
symmetric components (streamwise T x and wall- normal T y ) have significantly longer tails than a 
Gaussian function, and consequently exhibit higher than Gaussian values of kurtosis (cf. table [3]). 
The kurtosis is found to decrease with increasing particle size, consistent with the above filtering 
argument (also cf. figure [Tib) . 

The pdf of the fluctuations of the spanwise component of particle torque, T z , is clearly asym- 
metric with a pronounced negative skewness. Now, it is well established that the pdf of stream- 
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Figure 17: Normalised pdfs of the quantities for which statistical moments have been shown in 
figure |TH] The panels (a) and (c) show data from case F10 compared to data from the simple 
model with s + = 12; panels (b) and (d) show data from case F50 compared to data from the 
simple model with s + = 52. In Panels (a) and (b) the lines and symbols correspond to: dashed 
line, T^/cf,; dotted line, T^/erf,; Q, — -T^/ojt; —J'' z / a j r (please note the negative signs). Panels 
(c) and (d) show: dash-dotted line, T'/aj.; V, Tl/a^. The solid line corresponds to a Gaussian 
distribution. 



wise velocity fluctuations v! in smooth-w all channel flow is positively skewed close to the wall 
(|Kim et aU ' ll987t |j imenez fc Hovasl 120081) . In the limit of a wall-element with vanishing size, the 
pdf of T x is directly related to the pdf of the streamwise velocity fluctuations just above the wall. 
Since, as found above, the particle torque component T z behaves similarly as the smooth- wall force 
— T x (note the changed sign) , the observed negative skewness of the former is consistent with the 
positively skewed streamwise velocity pdf previously found in smooth-wall channel flow. 

It should be noted that the analogy drawn between shear forces acting upon a square element of 
a smooth-wall channel flow and the corresponding hydrodynamic torque components of spherical 
particles in the roughness layer can be expected to lose its appeal in the fully rough regime. In that 
case, which is outside the scope of the present study, pressure-induced forces will by far outweigh 
viscous forces. Although the torque around the particle center will still by definition be devoid of 
a pressure contribution (TfTj)) . other quantities of interest to the onset of particle motion, such as 
the torque around the line connecting the contact points with the downstream neighbor particles, 
might be dominated by the contribution from pressure forces. Therefore, the main utility of the 
proposed simple model is presumably limited to the regime of transitional roughness. 
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Figure 18: Non-normalized pdf of the lift coefficient defined as cl = F y /(pu 2 D 2 ir/4:), dashed line: 
case F10, dash-dotted line: case F50. 



3.5 Implications for the onset of sediment erosion 



Sediment erosion is often parametrised in terms of the non-dimensional Shields number 6 which 
is defined as follows: 

= 7 T " n , (14) 
(p p - p)gD 

where r w = pu 2 . is the wall shea r-stress, p n the density of the sedimen t part icles and g the value 
of the gravitational acceleration (Shi elds . 19361 Ivan Riinl . ll99aiGarciall2008h . If we suppose that 
erosion is initiated by lift forces alone, then we can characterise the onset of erosion by a balance 
between hydrodynamic lift force, F y , and buoyant weight of the particle, Fb = (p P — p)gD 3 ir/6, 
yielding the following expression for the critical Shields parameter: 



2 

3c L ' 



(15) 



where the lift coefficient is defined as cl = F y /{pvq.D tt/4). Please note, that in the above 
definition of a slightly different normalisation than in §3.31 is chosen. It can be seen that for 
this erosion scenario the critical value of the Shields number is inversely proportional to the lift 
coefficient at the onset of erosion. 

Figure [T5] shows the pdf of the lift coefficient Cl for the two present cases. In order to determine 
the smallest value of the Shields parameter for which sediment erosion can be initiated, the largest 
occurring value of needs to be considered in each case. Since the pdfs exhibit exponential tails, 
it is difficult to determine a precise upper bound of c/,. However, when a given (small) minimum 
probability of observation is fixed, it is clear that the larger spheres (case F50) will yield a larger 
maximum value of cl than the smaller spheres (F10). Consequently, a smaller critical Shields 
number is obtained for the spheres in case F50. Other modes of erosion (sliding parallel to the 
contact point with the downstream neighbour particles; rotation around the contact point) have 
also been analysed with similar conclusions (graphs of pdfs omitted) . 

Although exhibiting considerable scatter, experiments and field observations seem to indicate 
an increase with D + of the cri tical Shields number 9 C over the current range of particle diameters 
( van Riin . 1993 : Garcia . 2008). One has to be cautious, however, since results obtained in our 
idealised flow configuration are compared to experimental observations in flows involving a wide 
range of different irregular particle arrangements as well as varying particle shapes and size distri- 
butions. With this caveat in mind, we can take the different trend found in the present simulations 
(9 C decreases from D + = 10 to D + = 50) as compared to experiments (9 C increases from D + rj 10 
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to D + w 100) as an indication that extreme force- and torque-generating events recorded in fixed- 
particle configurations might not immediately yield a criterion which is sufficient to judge whether 
erosion will indeed occur as predicted when particles are freely mobile (under otherwise identical 
conditions). In particular, the question of the influence of the interaction between the incipient 
particle motion and the surrounding flow field as well as the related question of the necessary 
duration of force-/torque-generating flow events cannot be answered with certainty based upon 
fixed-particle data alone. In order to answer these queries, additional studies involving an analysis 
of high-fidelity data of the actual process of sediment erosion need to be carried out. 

4 Summary and conclusions 

Direct numerical simulation of open channel flow over a geometrically rough wall has been per- 
formed at a bulk Reynolds number of Reb w 2900. The wall consisted of a layer of spheres in 
a square arrangement touching a solid wall. Two particle diameters were considered: case F10 
with D+ = 10.7 (Re T = 188), and case F50 with D+ = 49.3 {Re T = 235). In case F10 the 
effect of roughness on the flow field statistics was small, and the limit of the hydraulically smooth 
flow regime is approached; in case F50 the roughness effect was stronger, and the flow is in the 
transitionally rough flow regime. 

The complexity of the time-averaged three-dimensional flow field within the roughness layer 
was discussed in detail. In both cases a recirculation forms downstream of the spheres that is 
connected over the entire spanwise direction, being more pronounced in the large sphere case. 
Three-dimensionality above the roughness layer is lost rapidly with wall-distance, yielding a time 
averaged flow field which is essentially one-dimensional beyond a distance of two particle diameters. 

The main result of this paper is the characterisation of the force and torque acting on a particle 
due to the turbulent flow. It was found that in the present cases the mean drag on a sphere is 
4% (case F10) and 15% (case F50) higher than the reference force Fr = pu^An, where An is the 
wall-normal projected area of the wall per particle. Given our definition of the friction velocity, 
these numbers reflect the fact that the drag force on the roughened bottom wall (below the fixed 
spheres) is small. In both cases a strong positive lift was obtained in agreement with previous 
experiments, exceeding values of 18% and 32%, respectively, of the corresponding drag. The values 
of the mean spanwise torque on a particle are comparable to —Frtr, where t\r is the distance 
from the particle centre to the position of the virtual wall, located at a distance of yo = 0-8D from 
the plane part of the solid wall. It was shown that in both cases the mean drag, lift and spanwise 
torque are to a large extent produced in a region of the particle surface which is located above the 
virtual wall [y > yo). The spatial distribution over the particle surface of the stresses that lead to 
time-averaged drag, lift and spanwise torque are found to be similar in shape in the two cases. 

We have observed that the intensity of particle force fluctuations (when normalised by Fr) is 
significantly larger in the large-sphere case. Conversely, when analysing the torque it is found that 
only the fluctuation intensity of the spanwise component is larger in the large-sphere case, whereas 
the two remaining components exhibit smaller fluctuation intensities when the sphere is larger. By 
means of a simplified model we were able to show that the torque fluctuations might be explained 
by the spheres acting as a filter with respect to the size of the flow scales which can effectively 
generate torque fluctuations. As a model we have considered the shear-forces and torque exerted 
by the flow on a square wall-element in a smooth-wall configuration. By systematically varying 
the linear dimension of the wall-element we are able to analyse the influence of the length scale. 
Here we find that the normalised fluctuation intensity of the streamwise and spanwise shear-forces 
monotonously decreases with the filter size, while the wall-normal torque experiences a maximum 
of normalised fluctuation intensity for intermediate filter sizes of approximately 70 wall units. By 
assuming that the largest part of the particle torque fluctuations is generated in a small area 
around the particle tops, the results from the simplified model carry over to the corresponding 
components of the particle torque. We obtain indeed a reasonable agreement between standard 
deviation and kurtosis of shear-forces and torque acting on a square wall-element on the one 
hand and respective particle torque components on the other hand. However, since we have only 
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considered two particle sizes, we cannot state with certainty that wall-normal particle torque 
fluctuations are indeed most intense at the above mentioned scale of 70 wall units. Similarly, 
based on the current data it is not possible to judge whether the model is capable of providing 
insight in the fully rough flow regime. These points should be clarified in future studies. 

Fluctuations of both force and torque were found to exhibit strongly non-Gaussian pdfs with 
particularly long tails. The deviation from a Gaussian distribution was significantly smaller in 
the large-sphere case, which was attributed to the smaller effect that highly intermittent small 
scales have on the larger particle surface area. Moreover, it was observed that the spanwise torque 
component has a marked negative skewness. In the light of the analogy with a wall-element in a 
smooth-wall configuration, this finding is consistent with the positive skewness of the streamwise 
velocity fluctuations near the wall in a smooth wall channel flow. 

Concerning the potential for sediment erosion, it was concluded from the present data that 
the largest recorded values of hydrodynamic lift lead to critical Shields numbers which are smaller 
in the large-sphere case as compared to the small-sphere case. The same trend was found when 
considering the forces projected onto the direction tangential to the downstream contact point 
between spheres in neighbouring positions as well as when evaluating the balance of angular 
moments around the contact point. Measurements in experiments with truly mobile particles 
seem to indicate the opposite trend (increasing critical Shields number with increasing particle 
size in the range of D + ss 10 to 100). However, these opposite trends do not necessarily imply a 
contradiction, since additional effects which might play a role in the dynamical process of erosion 
have not been addressed in the present study. 

Two important idealisations with respect to real- world sediment erosion have been made in the 
present work: the regularity of the geometrical arrangement and the immobility of the particles. 
Concerning the geometry, it is expected that different particle arrangements, size distributions 
and shapes will lead to a modification of the forces acting upon sediment particles. In particular, 
varying the protr usion of individual particles has be e n sho wn to have a significant effect on the 
onset of erosion ( Fenton fc Abbott . 1977 : Cameron . 20061 ). In this respect, the configuration 



studied in the present work can be considered as a case where mutual sheltering of particles is 
high due to their uniform diameter, spherical shape and regular arrangement. 

Concerning the immobility of the particles, we see several consequences arising from this ide- 
alisation which could potentially affect the implications for sediment erosion: (i) the modification 
of the flow field by the particles during the incipient motion; (ii) the determination of the tem- 
poral duration of force- and torque-generating flow events which is necessary in order to achieve 
irreversible onset of particle motion; (iii) the influence of collective mobility. In order to evaluate 
the importance of these mobility effects, additional data from configurations with truly eroding 
particles needs to be analysed. It might then still be possible to devise a refined erosion criterion 
which allows for the determination of erosion probabilities based upon data from fixed-particle 
configurations. 

Both of these additional aspects (the bed geometry and the particle mobility) should be ad- 
dressed in future studies. 
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of the manuscript. 



A The position of the virtual wall and the friction velocity 

Several common methods exist for the definition of an origin yo of the wall-normal coordinate 
when analysing turbulent flow statistics over rough walls. 

A priori definitions can be based on geometrical consideratio ns. Examples are the volume of 
the roughness elements divided by the area of the virtual wall (cf . ISchlichting . 19361 ) , which for the 
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present geometry leads to yo/D = 0.44 (0.56) in case F10 (F50), or the average of the maximum 
surface elevation, which leads to yo/D = 0.54 (0.65) in case F10 ( F50) . A poste riori methods 
emplo y the data from measurements or simulations to define yo- iThom (1971) and Jackson 
(Il98lh propose to define yo by the wall-normal position of the centroid of the drag profile on the 
roughness elements. In the present study such a definition would lead to values of yo/D — 0.88 
(0.84) in case F10 (F50). It should be noted that in case of a porous sediment layer, this definition 
is biased by the inter-porous flow. Most researchers, however , use methods which involve the 
adjustment of a logarithmic law to the mean velocity profile (jRaupach et especially 
for high Reynolds number flows. Based on these methods, several studies on turbulent flow over 
spherical roughness (for various Reynolds numbers, particle arrangements and flow ge ometries) 
can be found that provide the value of yp for a give n parti cle diameter ( cf. rev iews in iBavazit , 
Il983t iNezu fc Nakagawal Il993l p. 26; iDittrichl . Il998l p. 2 9; iDetert et all l2010al), including also 
studies that match well with the presen t flow conditions ( Nakagawa fc Nezul . Il977 ; Grass et "all . 
1991 ; Cameron . 20061 ; Isingh et ai , 12007 ). In these studies the virtual wall is positioned at yo/D in 



the interval 0.61 to 0.82. In som e studies the v i rtual wall is defined at the location of the roughness 
crests (i.e. yo/D 



1) such as in Manes et al 



(120071) . 

In the present work we have selected to fix the position of the virtual wall at a given level 
yo/D = 0.8 inside the range of values determined in relevant experiments. 

Turning now to the definition of the velocity scale u T , we will discuss three common approaches 
in the following. Again, a widely used method is to obtain u T by adjusting a logarithmic law to the 
mean velocity profile. Assuming the values k — 0.40 and yo/D — 0.80 for the Karman constant 
and the offset of the virtual wall, respectively, a fit over the range b05 v < (y — y ) < 0-5h yields 
Ur/Ubh = 0.062 (0.081) in case F10 (F50). However, it should be recalled that in the present low- 
Reynolds number flow the limited extent of the logarithmic region makes this approach relatively 
error-prone. 

Alternatively, the global momentum balance can be used in order to relate the driving force 
(either due to a pressure gradient or gravity) to the different contributions to the drag force 
generated at the fluid-solid interfaces. While the mean momentum balance is uniquely defined, it 
does not immediately provide a velocity scale. In some studies the velocity scale is defined from 
the volumetric for c e integrated from the virtu al wall-distance to the free surface (for example 
Nakagawa fc Nezu . 119771 IDetert et all 12010a ). i.e. in our notation = —(dpi/dx)h/p. This 
definition leads to u T /Ubh = 0.066 (0.081) in case F10 (F50). Ot her authors choose to integrate 
the driving force exclusively over the volume occupied by the fluid ( Manes et aZ.I . [2oTlh . leading to 



u^. A w = —(dpi/dx}Vf/p (where A w is the area of the wall-parallel cross-section of the considered 
control volume and Vf the corresponding volume occupied by the fluid). This latter definition 
yields u T /U bh = 0.066 (0.081) in case F10 (F50). 

Finally, let us consider definitions based on the total shear stress profile. In smooth-wall flow, 
the total shear stress r tot is linear with wall-distance and the appropriate velocity scale is given 
by = T to t(0)/p. In rough- wall flow, T to t in general deviates from a linear relation below the 
roughness crests which prevents the use of a similar definition, e.g. based upon Tt t{y = yo)- 
Instead, some researchers propose to determine the velocity scale independently of the position of 
the virtual wall by us ing the total shear stress at the roughness crests, i.e. v? T — T tot (y — D)/p 



( Pokraiac et ai , 2006). This definition leads to values of u T /Ubh = 0.066 (0.080) in case F10 



(F50). Note that this latter definition makes a direct comparison of different data sets difficult, 
since the total shear stress profiles T to t/(pu^.) represented as a function of (y — yo)/h will in general 
not collapse. Alternatively, u T can be computed from the total shear stress extrapolated from the 
region where it varies linearly (i.e. above the roughness crests) down to the position of the virtual 
wall, yielding the defining relation 



2 I V ' 

r to t = pu T 1 



(16) 



valid for y > D. This definition leads to u T /U bh = 0.066 (0.082) in case F10 (F50). Incidentally it 
can be deduced from the global momentum balance that our definition implies u^. = ~(dpi/dx)h/ p, 
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i.e. it tu rns out that the definition of u T throu gh (TTBl) is equiv alent to the above mentioned definition 
used by iNakagawa &; Nezu (1977) and Detert et al. ( 2010a ) based upon an integral of the driving 
force. 



B Details on averaging procedures 



In the present paper two definitions have been used for averaging a discrete flow field, <j>(i,j,k), 
in x-z planes: 



1 



N x N z 
1 



%=\ k=i 



22'^2<l>(i,j,k)m(i,j,k) , 



(17) 
(18) 



where N x and N z are the number of grid points in x and z directions, respectively, and m is a field 
that works as a mask for computing the averages. If a given point lies within the fluid domain, 
at that point m = 1, otherwise m = 0. N m (j) is the sum of m over a wall parallel plane at wall- 
distance yj, i.e. N m (j) — Yli=i Hfc=i m (hj> ty- N m (j) equals N X N Z above the roughness layer 
such that both expressions (fl"8|) and fI7| are equal to each other away from the roughness elements. 
Within the roughness layer N m (j) and thus the number of samples for each wall parallel plane 
decreases. In the context of the double-averaging methodology these t wo quantities are ge nerally 
referred to as superficial and intrinsic spatial average, respectively (cf. lNikora et al . 2007 ). 

Note that the zero velocity condition is forced only at the surface of the particles due to reasons 
of efficiency ( Uhlmannl . 2005a). This l e ads to fictitious non-zero velocities at the grid points that 
lie within the particles. Fadlun et al. ( 2000h demonstrated that t he extern a l flow i s essentially 
unchanged by this procedure which has been confirmed later by luhi mannl (|2005al ). Since the 
internal fictitious flow affects the value of ((j)} xz (j) (according to [T7]) in the roughness layer, we 
present averages computed according to (fT8|) where the flow within the roughness layer is discussed 
(i.e. figure [5] and figure [8]). When focusing upon the flow above the roughness layer (i.e. in figure [2] 
and figure BJ, we choose to present data computed according to (|17|) . because the number of 
available samples is larger, as explained in § O 



C Consistency of runtime and a posteriori statistics 

For the flow field statistics presented in this paper, two different sets of data have been used. 
The first set of flow field statistics was collected during the runtime of the simulation employing 
equation (fT7|) . This leads to a number of the order of 10 11 samples per wall-normal grid point 
in case F10 and F50, collected over the entire observation interval. The second set of data was 
obtained from analysing stored snapshots of the flow field of which 90 were used in each case. 
The latter set has been used to compute some additional statistical quantities not stored during 
runtime. Since it provides a smaller number of samples (roughly a factor of 10 3 less), we will in 
the following check its consistency with the more complete set accumulated at runtime. 

Figure [T5] shows for each case the second order moments of the velocity fluctuations obtained 
at run-time in comparison to the same quantities obtained from the snapshots of the simulations 
applying the averaging operator as defined in (fl7|). The differences between the two data sets are 
small, measuring less than 0.06u T (0.02u T ) in case F10 (F50). Incidentally, it can be seen from 
the figures that the discrepancy is largest near the open surface. We can therefore conclude, that 
the data set provided from the 90 stored snapshots is sufficient for the purpose of computing the 
quantities shown in figure [5l figure |6l figure [7] and figure [8j 
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Figure 19: Comparison of velocity fluctuations normalised by u T obtained from run-time (solid 
line) and from snapshots (dashed line) as a function of y/H. Curves from left to right are the 
components in wall- normal (v rms /u T ), spanwise (w rms /u T ) and streamwise {u rms /u T ) direction. 
Panel (a) shows case F10, panel (b) shows case F50. 
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